home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zheev.f < prev    next >
Text File  |  1997-06-25  |  7KB  |  206 lines
  1.       SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
  2.      $                  INFO )
  3. *
  4. *  -- LAPACK driver routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          JOBZ, UPLO
  11.       INTEGER            INFO, LDA, LWORK, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       DOUBLE PRECISION   RWORK( * ), W( * )
  15.       COMPLEX*16         A( LDA, * ), WORK( * )
  16. *     ..
  17. *
  18. *  Purpose
  19. *  =======
  20. *
  21. *  ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
  22. *  complex Hermitian matrix A.
  23. *
  24. *  Arguments
  25. *  =========
  26. *
  27. *  JOBZ    (input) CHARACTER*1
  28. *          = 'N':  Compute eigenvalues only;
  29. *          = 'V':  Compute eigenvalues and eigenvectors.
  30. *
  31. *  UPLO    (input) CHARACTER*1
  32. *          = 'U':  Upper triangle of A is stored;
  33. *          = 'L':  Lower triangle of A is stored.
  34. *
  35. *  N       (input) INTEGER
  36. *          The order of the matrix A.  N >= 0.
  37. *
  38. *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
  39. *          On entry, the Hermitian matrix A.  If UPLO = 'U', the
  40. *          leading N-by-N upper triangular part of A contains the
  41. *          upper triangular part of the matrix A.  If UPLO = 'L',
  42. *          the leading N-by-N lower triangular part of A contains
  43. *          the lower triangular part of the matrix A.
  44. *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  45. *          orthonormal eigenvectors of the matrix A.
  46. *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
  47. *          or the upper triangle (if UPLO='U') of A, including the
  48. *          diagonal, is destroyed.
  49. *
  50. *  LDA     (input) INTEGER
  51. *          The leading dimension of the array A.  LDA >= max(1,N).
  52. *
  53. *  W       (output) DOUBLE PRECISION array, dimension (N)
  54. *          If INFO = 0, the eigenvalues in ascending order.
  55. *
  56. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  57. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  58. *
  59. *  LWORK   (input) INTEGER
  60. *          The length of the array WORK.  LWORK >= max(1,2*N-1).
  61. *          For optimal efficiency, LWORK >= (NB+1)*N,
  62. *          where NB is the blocksize for ZHETRD returned by ILAENV.
  63. *
  64. *  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
  65. *
  66. *  INFO    (output) INTEGER
  67. *          = 0:  successful exit
  68. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  69. *          > 0:  if INFO = i, the algorithm failed to converge; i
  70. *                off-diagonal elements of an intermediate tridiagonal
  71. *                form did not converge to zero.
  72. *
  73. *  =====================================================================
  74. *
  75. *     .. Parameters ..
  76.       DOUBLE PRECISION   ZERO, ONE
  77.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  78.       COMPLEX*16         CONE
  79.       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
  80. *     ..
  81. *     .. Local Scalars ..
  82.       LOGICAL            LOWER, WANTZ
  83.       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
  84.      $                   LLWORK, LOPT
  85.       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
  86.      $                   SMLNUM
  87. *     ..
  88. *     .. External Functions ..
  89.       LOGICAL            LSAME
  90.       DOUBLE PRECISION   DLAMCH, ZLANHE
  91.       EXTERNAL           LSAME, DLAMCH, ZLANHE
  92. *     ..
  93. *     .. External Subroutines ..
  94.       EXTERNAL           DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR,
  95.      $                   ZUNGTR
  96. *     ..
  97. *     .. Intrinsic Functions ..
  98.       INTRINSIC          MAX, SQRT
  99. *     ..
  100. *     .. Executable Statements ..
  101. *
  102. *     Test the input parameters.
  103. *
  104.       WANTZ = LSAME( JOBZ, 'V' )
  105.       LOWER = LSAME( UPLO, 'L' )
  106. *
  107.       INFO = 0
  108.       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  109.          INFO = -1
  110.       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
  111.          INFO = -2
  112.       ELSE IF( N.LT.0 ) THEN
  113.          INFO = -3
  114.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  115.          INFO = -5
  116.       ELSE IF( LWORK.LT.MAX( 1, 2*N-1 ) ) THEN
  117.          INFO = -8
  118.       END IF
  119. *
  120.       IF( INFO.NE.0 ) THEN
  121.          CALL XERBLA( 'ZHEEV ', -INFO )
  122.          RETURN
  123.       END IF
  124. *
  125. *     Quick return if possible
  126. *
  127.       IF( N.EQ.0 ) THEN
  128.          WORK( 1 ) = 1
  129.          RETURN
  130.       END IF
  131. *
  132.       IF( N.EQ.1 ) THEN
  133.          W( 1 ) = A( 1, 1 )
  134.          WORK( 1 ) = 3
  135.          IF( WANTZ )
  136.      $      A( 1, 1 ) = CONE
  137.          RETURN
  138.       END IF
  139. *
  140. *     Get machine constants.
  141. *
  142.       SAFMIN = DLAMCH( 'Safe minimum' )
  143.       EPS = DLAMCH( 'Precision' )
  144.       SMLNUM = SAFMIN / EPS
  145.       BIGNUM = ONE / SMLNUM
  146.       RMIN = SQRT( SMLNUM )
  147.       RMAX = SQRT( BIGNUM )
  148. *
  149. *     Scale matrix to allowable range, if necessary.
  150. *
  151.       ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
  152.       ISCALE = 0
  153.       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
  154.          ISCALE = 1
  155.          SIGMA = RMIN / ANRM
  156.       ELSE IF( ANRM.GT.RMAX ) THEN
  157.          ISCALE = 1
  158.          SIGMA = RMAX / ANRM
  159.       END IF
  160.       IF( ISCALE.EQ.1 )
  161.      $   CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
  162. *
  163. *     Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
  164. *
  165.       INDE = 1
  166.       INDTAU = 1
  167.       INDWRK = INDTAU + N
  168.       LLWORK = LWORK - INDWRK + 1
  169.       CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
  170.      $             WORK( INDWRK ), LLWORK, IINFO )
  171.       LOPT = N + WORK( INDWRK )
  172. *
  173. *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
  174. *     ZUNGTR to generate the unitary matrix, then call ZSTEQR.
  175. *
  176.       IF( .NOT.WANTZ ) THEN
  177.          CALL DSTERF( N, W, RWORK( INDE ), INFO )
  178.       ELSE
  179.          CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
  180.      $                LLWORK, IINFO )
  181.          INDWRK = INDE + N
  182.          CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA,
  183.      $                RWORK( INDWRK ), INFO )
  184.       END IF
  185. *
  186. *     If matrix was scaled, then rescale eigenvalues appropriately.
  187. *
  188.       IF( ISCALE.EQ.1 ) THEN
  189.          IF( INFO.EQ.0 ) THEN
  190.             IMAX = N
  191.          ELSE
  192.             IMAX = INFO - 1
  193.          END IF
  194.          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
  195.       END IF
  196. *
  197. *     Set WORK(1) to optimal complex workspace size.
  198. *
  199.       WORK( 1 ) = MAX( 2*N-1, LOPT )
  200. *
  201.       RETURN
  202. *
  203. *     End of ZHEEV
  204. *
  205.       END
  206.